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Abstract 

We  design  fast  protocols  to  separate  or  recombine  two  ions  in  a  segmented  Paul  trap.  By  inverse 
engineering  the  time  evolution  of  the  trapping  potential  composed  of  a  harmonic  and  a  quartic  term,  it  is 
possible  to  perform  these  processes  in  a  few  microseconds  without  final  excitation.  These  times  are  much 
shorter  than  the  ones  reported  so  far  experimentally.  The  design  is  based  on  dynamical  invariants  and 
dynamical  normal  modes.  Anharmonicities  beyond  the  harmonic  approximation  at  potential  minima 
are  taken  into  account  perturbatively.  The  stability  versus  an  unknown  potential  bias  is  also  studied. 


1.  Introduction 


Trapped  cold  ions  provide  a  leading  platform  to  implement  quantum  information  processing.  Separating  ion  chains 
is  in  the  toolkit  of  basic  operations  required.  (Merging  chains  is  the  corresponding  reverse  operation  so  we  shall  only 
refer  to  separation  hereafter.)  It  has  been  used  to  implement  two-qubit  quantum  gates  [1];  also  to  purify  entangled 
states  [2, 3],  or  teleport  material  qubits  [4].  Moreover,  an  architecture  for  processing  information  scalable  to  many 
ions  could  be  developed  based  on  shuttling,  separating  and  merging  ion  crystals  in  multisegmented  traps  [5]. 

Ion-chain  separation  is  known  to  be  a  difficult  operation  [6].  Experiments  have  progressed  towards  lower 
final  excitations  and  shorter  times  but  much  room  for  improvement  still  remains  [7-9].  Problems  identified 
include  anomalous  heating,  so  devising  short-time  protocols  via  shortcuts  to  adiabaticity  (STA)  techniques  was 
proposed  as  a  way-out  worth  exploring  [10].  STA  methods  intend  to  speed  up  different  adiabatic  operations 
[11,12]  without  inducing  final  excitations.  An  example  of  an  elementary  (fast  quasi- adiabatic)  STA  approach 
[11]  was  already  applied  for  fast  chain  splitting  in  [7].  Here,  we  design,  using  a  more  general  and  efficient  STA 
approach  based  on  dynamical  normal  modes  (NM)  [13,14],  protocols  to  effectively  separate  two  equal  ions, 
initially  in  a  common  electrostatic  linear  harmonic  trap,  into  a  final  configuration  where  each  ion  is  in  a  different 
well.  The  motion  is  assumed  to  be  effectively  one-dimensional  due  to  tight  radial  confinement.  The  external 
potential  for  an  ion  at  q  is  approximated  as 

^ext  =  a(t)q2  +  P(t)q\  (1) 

which  is  experimentally  realizable  with  state-of-art  segmented  Paul  traps  [6, 15]. 

Using  dynamical  NM  [13, 14]  a  Hamiltonian  will  be  set  which  is  separable  in  a  harmonic  approximation 
around  potential  minima.  By  means  of  Lewis-Riesenfeld  invariants  [16]  we  shall  design  first  the  approximate 
dynamics  of  an  unexcited  splitting,  taking  into  account  anharmonicities  in  a  perturbative  manner,  and  from  that 
inversely  find  the  corresponding  protocol,  i.e.  the  a(t)  and  /3(t)  functions. 

The  Hamiltonian  of  the  system  of  two  ions  of  mass  m  and  charge  e  is,  in  the  laboratory  frame 


2m  2m 


V  —  a(t)(^q^  +  g22)  +  (3(t)(q^  +  q^  H - - — >  (2) 

%  ~  *?2 
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Figure  1 .  Scheme  of  the  separation  process.  At  t  =  0  (left),  both  ions  are  trapped  within  the  same  external  harmonic  potential.  At  final 
time  tf  (right),  the  negative  harmonic  term,  and  a  quartic  term  build  a  double  well  external  potential.  The  aim  of  the  process  is  to  set 
each  of  the  ions  in  a  different  well  without  any  residual  excitation. 


wherep!,p2  are  the  momentum  operators  for  both  ions,  qu  q2  their  position  operators,  and  Q  —  e0  being 
the  vacuum  permittivity.  We  use  on  purpose  a  c-number  notation  since  we  shall  also  consider  classical 
simulations.  The  context  will  make  clear  if  c-numbers  or  ^-numbers  are  required.  We  suppose  that,  due  to  the 
strong  Coulomb  repulsion,  qx  >  qv  By  minimizing  the  potential  part  of  the  Hamiltonian  V,  we  find  for  the 
equilibrium  distance  between  the  two  ions,  d(t),  the  quintic  equation  [6] 

(3(t)d5(t)  +  2 a(t)d\t)  -  2 Cc  =  0,  (3) 

which  will  be  quite  useful  for  inverse  engineering  the  ion-chain  splitting,  even  without  an  explicit  solution  for  d 
( t ).  At  t  =  0  a  single  external  well  is  assumed,  (3(0)  =  0  and  a  (0)  >  0,  whereas  in  the  final  double-well 
configuration >  0,  a(t{)  <  0.  At  some  intermediate  time  tj  the  potential  becomes  purely  quartic 
( a  (ti )  =  0).  Our  aim  is  to  design  the  functions  a(t)  and  / 3(t )  so  that  each  of  the  ions  ends  up  in  a  different 
external  well  as  shown  in  figure  1,  in  times  as  short  as  possible,  and  without  any  final  excitation. 


2.  Dynamical  Normal  Modes 


To  define  dynamical  NM  coordinates,  we  calculate  first  at  equilibrium  (the  point  {q!°\  q 2(0) }  in  configuration 


space  where  the  potential  is  a  minimum,  dV /dq1  =  dV /dq2  =  0)  the  matrix  Vy  =  — 


1  ’  12 
1  d2V 


I  ea  ’ 


m  dqt  dq-  eci 

The  equilibrium  positions  are  q^y>  =  and  the  matrix  takes  the  form 

\ 


1 


2  J2  2 

2  Cc  2  Cc 

2a  +  12/3—  H - y  - 7- 


m 


4  d 3  d? 

-2  Q  ^  ^ad2  2Cc 

— 2a  +  12/3—  H - 

d 3  4  d3 


(4) 


The  eigenvalues  are 


A_  =  —  (2  a  +  3  /3d2\ 
m  v  J 


A+  —  —  f  2a  +  3 (3d2  -f — -2 

m  \  d J 


1  r 


(5) 


which  define  the  NM  frequencies  as  Q±  =  y[\±  corresponding  to  center-of-mass  (— )  and  relative  (stretch) 
motions  (+).  These  relations,  with  equation  (3)  written  as 


P(t)  = 


2  Cc  2  a(t) 


d5(t)  d2(t) 

allow  us  to  write  a(t)  and  d(t)  as  functions  of  the  NM  frequencies 

a(t)  =  — 


(6) 


(7) 


d  (t)  = 


4  C 


(8) 


Substituting  these  expressions  into  equation  (6),  /3(t)  may  also  be  written  in  terms  of  NM  frequencies. 
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The  normalized  eigenvectors  are 


JrCl 


V+: 


which  we  denote  as  v± 

1  a  bo  ra  t  o  ry  coo  rd  i  n  a  t  es  as 


(:;)• 


MW 


(9) 


The  (mass  weighted)  dynamical  NM  coordinates  are  defined  in  terms  of  the 

q±  =  fl±Vm  (<?!  -  ^(0))  +  b±~Jm(q2  -  q 2(0)). 


(10) 


The  unitary  transformation-of-coordinates  matrix  is 

u  =  (q+,  q_| qv  q2)  =  -  <fi(q+,  q_)]<fi[^2  “  «2(<l+>  I-)]'  (11) 

The  Hamiltonian  in  the  dynamical  equation  for  |  vp')  =  U\ip),  where  |  ip>)  is  the  lab-frame  time-dependent  wave 
function  (in  qu  q2  coordinates)  evolving  with  H ,  is  given  by 


H'  =  UHU f  -  iHu(dtW) 

P+  1  (\1  2  d 

~T  + 2  +q++ vT^P+ 

+  —  +  q2  5  (12) 
2  2 


plus  qubic  and  higher  order  terms  in  the  potential  that  we  neglect  by  now  (they  will  be  considered  in  section  4 
below).  Similarly  to  [13, 14],  we  apply  a  further  unitary  transformation  U  =  e_K^q+/ to  write  down  an 
effective  Hamiltonian  for  |  ipff)  =  U\  'ip1)  with  the  form  of  two  independent  harmonic  oscillators  in  NM  space, 
H"  =  UH'W  -  i m(dtU\ 


H"  =  —  + 

2  2 


Jmd 

q+  +  VfiT 


—  +  -SlW  =  H l  +  H". 

2  2“ 


(13) 


These  oscillators  have  dynamical  invariants  of  the  form  [16] 


[p±(p±  -  *±)  -  At(q±  -  *±)J 


+ 


— fl 


o± 


/  \2 
q^-x± 

P± 


(14) 


where  the  auxiliary  functions  p±  and  x+  satisfy 


At  +  $d±P±  — 

x+  +  Q2+x+  =  - 


(15) 

(16) 


with  Q0±  =  £2±(0),  and,  due  to  symmetry,  %_  =  0. 

The  physical  meaning  of  the  auxiliary  functions  maybe  grasped  from  the  solutions  of  the  time- dependent 
Schrodinger  equations  (for  each  NM  Hamiltonian  H±  in  equation  (13))  proportional  to  the  invariant 
eigenvectors  [22].  They  form  a  complete  basis  for  the  space  spanned  by  each  Hamiltonian  and  take  the  form 


Kd l(0>  = 


en  L 


^  +  (x±P±-x±p±) 


$n(o±) 


PL/2 


(17) 


where  cr±  =  — — — ,  and  (o±)  are  the  eigenfunctions  of  the  static  harmonic  oscillator  at  time  t—  0.  Thus  p±  are 

P i 

scaling  factors  proportional  to  the  state  width’  in  NM  coordinates,  whereas  the  x±  are  the  dynamical-mode 
centers  in  the  space  of  NM  coordinates.  Within  the  harmonic  approximation  there  are  dynamical  states  of  the 
factorized  form  |  pj"  (£) )  =  |  ip+  ( t ) )  |  ip'!_  ( t ) )  for  the  ion  chain  dynamics,  where  the  NM  wave  functions  |  ip±  ( t ) ) 
evolve  independently  with  H± .  They  may  be  written  as  combinations  of  the  form  |  ip>±(t))  =  Enc„±\€±m 
with  constant  amplitudes  Cn±.  The  average  energies  of  the  nth  basis  states  for  the  two  NM  are 


3 


IOP  Publishing 


New  J.  Phys.  17(2015)093031 


K±  =  «±l  ff±K±>, 


(2n  H-  l)/i 


-cn+ 


4^o- 
(In  +  l)/i 


e?_  =  ^  +  0 


4f2, 


o+ 


P+  +  &+P+  + 


P2. 

n20+' 


1  .2 

H — Xi 
2 


1  o2  I  ..  VOTrf 


2fi+r+  vin* 
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(18) 


3.  Design  of  the  control  parameters 


Once  the  Hamiltonian  and  Lewis-Riesenfeld  invariants  are  defined,  we  proceed  to  apply  the  invariant  based 
inverse  engineering  technique  and  design  ST  A.  The  results  for  the  simple  harmonic  oscillator  in  [1 1]  serve  as  a 
reference  but  have  to  be  extended  here  since  the  two  modes  are  not  really  independent  from  the  perspective  of 
the  inverse  problem.  This  is  because  a  unique  protocol,  i.e.,  a  single  set  of  a(t)  and  (3(f)  functions  has  to  be 
designed. 

We  first  set  the  initial  and  target  values  for  the  control  parameters  a(t)  and  / 3(t ).  At  time  t—  0,  the  external 
trap — for  a  single  ion — is  purely  harmonic,  with  (angular)  frequency  cj0.  From  equation  (5),  we  find  that 

FL(0)  =  ujq  and  f2+(0)  =  uq.  The  equilibrium  distance  is  d  (0)  =  •  For  the  final  time,  we  set  a  tenfold 

expansion  of  the  equilibrium  distance,  d(tf)  =  10d(0),  and  Q-(t{)  =  cjq.  This  also  implies 

F2+(tf )  =  V  1.002  u  o  ~  Q-(t{ ),  i.e.,  the  final  frequencies  of  both  NM  are  essentially  equal,  the  Coulomb 

interaction  is  negligible,  and  the  ions  can  be  considered  to  oscillate  in  independent  traps. 

The  inverse  problem  is  somewhat  similar  to  the  expansion  of  a  trap  with  two  equal  ions  in  [14],  but 
complicated  by  the  richer  structure  of  the  external  potential.  The  Hamiltonian  (13)  and  the  invariant  (14)  must 
commute  at  both  boundary  times  [H(tb),  /(%)]  =  0, 

th  =  0,  tf,  (19) 


to  drive  initial  levels  into  final  levels  via  a  one-to-one  mapping.  This  is  achieved  by  applying  appropriate 
boundary  conditions  (BCs)  to  the  auxiliary  functions  p± ,  x±  and  their  derivatives 

p±(0)  =  1,  P±(tf)  =  7±>  (20) 

P±(fb)  =  P±(fb)  =  °’  (2V> 

x+(fb)  =  *+(fb)  =  x+O  =  0,  (22) 

where  7±  =  J  777^7 .  Let  us  recall  that  x_  =  0  for  all  times  so  this  parameter  does  not  have  to  be  considered 
further. 

Inserting  the  BC  for  x+fo)  and  x+fo)  in  equation  (16)  we  find  that  d(%)  =  0.  Additionally,  d  (t\>)  =  0  is  to 
be  imposed  so  that  IA(%)  =  1.  According  to  equation  (8),  d  (%)  =  d(%)  =  0  by  imposing 
f2±(tb)  =  f^±(*b)  =  0.  With  d(t b)  =  d(t b)  =  0  the  Hamiltonians  and  wave  functions  coincide  at  the 
boundaries,  H^tb)  =  H"(%),  \fr  (%))  =  l^^b))*  which  simplifies  the  calculation  of  the  excitation  energy. 
From  equation  (15),  the  NM  frequencies  maybe  written  as 


Cl±  = 


At 
p±  ' 


(23) 


Thus  the  BC  f2±(tb )  =  &±(tb )  =  0  are  satisfied  by  imposing  on  the  auxiliary  functions  the  additional  BC 

]9±(fb)  =  'P±(fb)  =  0.  (24) 

We  may  now  design  ansatzes  for  the  auxiliary  functions  p±  that  satisfy  the  ten  BC  in  equations  (20),  (21)  and 
(24),  plus  the  BC  for  x+(%)  and  x+(%)  in  equation  (22)  (since  d  (t^)  =  0,  x+(%)  =  0  is  then  automatically 
satisfied,  see  equation  (16)).  Finally,  from  the  NM  frequencies  given  by  equation  (23)  we  can  inverse  engineer  the 
control  parameters  a  (t)  and  (3  (t)  from  equations  (6),  (7)  and  (8) . 

A  simple  choice  for  p_(t)  is  a  polynomial  ansatz  of  9th  order  p_  =  yT_o  bjs\  where  s  =  t/tf.  Substituting 
this  form  in  the  ten  BC  in  equations  (20),  (21)  and  (24),  we  finally  get 
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tf  (a/s)  tf  (a/s)  tf  (a/s) 

Figure  3.  (a)  Final  excitation  energy  versus  final  time  using  the  inverse  engineering  design  of  section  3  (solid  blue),  and  the  design  that 
takes  into  account  anharmonicities  in  section  4  (dashed  red),  (b)  Values  of  the  free  parameters  a10  (solid  blue)  and  ax  1  (dashed  red)  that 
minimize  the  excitation  energy  for  the  11th  order  polynomial  (Al).  (c)  Parameters  Ci0  (solid  blue),  cn  (dashed  red)  and  c12  (dashed-dotted 
green)  that  minimize  the  excitation  energy  for  the  1 2th  order  polynomial  (A2).  Two  9Be+  ions  where  split,  with  uj0 / (2tt)  =  2  MHz. 


p_  =  126(7 -  -  l)s5  -  420(7-  -  l)s6  +  540(7-  -  l)s7 

-  315(7-  -  l)s8  +  70(7-  -  l)s9  +  1.  (25) 

For  p+  we  will  use  an  11th  order  polynomial  p+  =  ansn  to  satisfy  as  well  x+(t\>)  =  x+(%)  =  O.The 

parameters  a0_ 9  are  fixed  so  that  the  ten  BC  for  p+  are  fulfilled  (see  the  appendix),  whereas  a10,  an  are  left  free, 
and  will  be  numerically  determined  by  a  shooting  program  [  17]  ('fminsearch’  in  MATLAB  which  uses  the 
Nelder-Mead  simplex  method  for  optimization),  so  that  the  remaining  BC  for  x+(%)  and  x+(%)  are  also 
satisfied.  Specifically,  for  each  pair  { a10,  au  }>  0±(t)  and  d(t)  are  determined  from  equations  (8)  and  (15),  to 
solve  equation  (16)  for  x+(t)  with  initial  conditions  x+(0)  =  x+(0)  =  0.  The  free  constants  are  changed  until 
x+(tf )  =  0  and  x+(tf )  =  0  are  satisfied.  Numerically  a  convenient  way  to  find  the  solution  is  to  minimize  the 
energy  E"+  (tf )  in  equation  ( 1 8). 

Figures  2(a)  and  (b)  depict  the  control  parameters  a  (t)  and  /?  (t)  found  with  this  method,  using 
equations  (6)  and  (7),  for  some  value  of  tf  and  ujq,  see  the  caption,  while  figure  2(c)  represents  the  equilibrium 
distance  between  ions  as  a  function  of  time  (8),  and  figure  2(d)  the  NM  frequencies.  In  figure  3(a)  the  excitation 
energy  is  shown,  versus  final  time,  for  optimized  parameters  given  in  figure  3(b).  The  initial  state  is  the  ground 
state  of  the  two  ions.  It  is  calculated  by  propagating  an  initial  guess  of  the  wave  function  in  imaginary  time  until  it 
relaxes  to  the  lowest  eigenfunction  [18].  The  excitation  energy  is  Eex  =  E(tf)  —  E0(tf),  where  E  ( tf )  is  the  final 
energy,  calculated  in  the  lab  frame,  and  E0  (tf )  is  the  final  ground-state  energy.  The  wave  function  evolution  is 
calculated  using  the  'split-operator  method’  with  the  full  Hamiltonian  (2).  If  the  harmonic  approximation  were 
exact,  there  would  not  be  any  excitation  with  this  STA  method,  E  ( tf )  =  Eq+  ( tf )  +  Eq_  ( tf )  =  E0  (tf ),  see 
equation  (18).  The  actual  result  is  perturbed  by  the  anharmonicities  and  NM  couplings.  The  final  ground  state  is 
also  calculated  with  an  'imaginary- time  evolution’.  The  corresponding  final  ground  state  energy  is  essentially 
two  times  a  harmonic  oscillator  ground  state  energy  plus  the  (negligible)  Coulomb  repulsion  at  distance  d(tf ). 
For  the  final  times  of  all  the  examples,  as  it  was  noted  in  previous  works  [10, 13, 14, 19],  classical  simulations 
(solving  Hamilton’s  equations  from  the  equilibrium  configuration  instead  of  Schrodinger’s  equation)  give 
indistinguishable  results  in  the  scale  of  figure  3(a). 

The  excitation  energy  in  figure  3(a)  (solid  line)  increases  at  short  times  since  the  harmonic  NM 
approximation  fails  [13,14].  However,  it  goes  down  rapidly  below  one  excitation  quantum  at  times  which  are 
still  rather  small  compared  to  experimental  values  used  so  far  [7, 8].  In  the  following  section  we  shall  apply  a 
perturbative  technique  to  minimize  the  excitation  further. 
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Table  1.  Maximum  values  of  /3,  and  critical  times 
(final  times  at  which  excitations  below  0. 1  quanta 
are  reached)  for  different  values  of  uj0.  The  calcu¬ 
lations  were  performed  with  the  11th  order  poly¬ 
nomial  for  p+. 


u0  (MHz) 

/WlO“3Nm-3) 

tcrit  (MS) 

3 

44.2 

2.9 

2 

11.4 

4.4 

1.2 

2.082 

7.4 

0.8 

0.539 

11.2 

4.  Beyond  the  harmonic  approximation 

An  improvement  of  the  protocol  is  to  consider  the  perturbation  of  the  higher  order  terms  neglected  in  the 
Hamiltonian  (13).  These  ‘anharmonicities’  [20]  are  cubic  and  higher  order  terms  in  the  Taylor  expansion  of  the 
Coulomb  term  Cc/(ql  —  q2 ), 

6V  =  Y°°JV® 

In  NM  coordinates  the  terms  take  the  simple  form 

6VV  =  (_1);+1^[ JTq  V;  (27) 

\  \  m  ) 

which  maybe  regarded  as  a  perturbation  to  be  added  to  H"  in  equation  (13).  (The  perturbation  does  not  couple 
the  center-of-mass  and  relative  subspaces.)  To  first  order,  the  excess  energy  due  to  these  perturbative  terms  at 
final  time  is  given  by 

toft.  =  (V4V(ff)l^°'>l€+(rf)}>  (28) 

where  the  |^+)  are  the  unperturbed  states  in  equation  (17).  Inverse  engineering  the  splitting  process  may  now 
be  carried  out  by  considering  a  12th  order  polynomial  for  p+  (see  (A2)),  with  three  free  parameters  so  as  to  fix  the 
BC  for  x+  and  also  minimize  the  excitation  energy.  In  practice  we  use  MATLAB’s  TminsearclT  function  for  the 
shooting  to  minimize  E0 +(tf )  +  8E^_  as  no  significant  improvement  occurs  by  including  higher  order  terms. 
Figure  3(a)  compares  the  performance  of  such  a  protocol  with  the  simpler  one  with  the  1 1th- order  polynomial 
(Al).  Figure  3(c)  gives  the  values  of  optimized  parameters  at  different  final  times. 

5.  Discussion 

A  large  quartic  potential  is  desirable  to  control  the  excitations  produced  at  the  point  where  the  harmonic  term 
changes  its  sign  [10].  At  this  point,  the  harmonic  potential  switches  from  confining  to  repulsive,  which  reduces 
the  control  of  the  system  and  potentially  increases  diabaticities  and  heating.  In  the  inverse  approach  proposed 
here  there  is  no  special  design  of  the  protocol  at  this  point,  but  the  method  naturally  seeks  high  quartic 
confinements  there.  In  figure  2(b)  (3  reaches  its  maximum  value  right  at  the  time  where  a  changes  sign  (see 
figure  2(a)).  However,  the  maximum  value  that  (3  can  reach  will  typically  be  limited  in  a  Paul  trap  [6]. 

In  table  1  we  summarize  the  different  maximal  values  of  (3  and  critical  times  (final  times  at  which  excitations 
below  0. 1  quanta  are  reached)  for  different  values  of  using  the  11th  order  polynomial  (Al)  for  p+. 

The  maximum  (3  decreases  with  tf,  such  that  the  shortest  possible  tf  at  a  given  maximum  tolerable  excitation 
energy  is  limited  by  the  achievable  / 3 .  The  trap  used  in  [8]  yields  a  maximum  f3  of  about  1 0-4  N  m-3,  at  =L  1 0  V 
steering  range.  In  a  recent  experiment  reported  in  [23],  the  value  used  was  (3  «  5  x  10-3  N  m-3.  The  numbers 
reported  in  the  table  are  thus  within  reach,  as  the  / 3  coefficients  scale  with  the  inverse  4th  power  of  the  overall  trap 
dimension,  and  technological  improvements  on  arbitrary  waveform  generators  may  allow  for  operation  at  an 
increased  voltage  range. 

Another  potential  limitation  the  method  could  encounter  in  the  laboratory  is  due  to  biases  (a  linear  slope)  in 
the  trapping  potential,  Vext  =  aq2  +  (3qA  +  \q,  with  A  constant  and  unknown  [10].  Figure  4  represents  the 
excitation  energy  versus  the  energy  difference  between  the  two  final  minima  of  the  external  potential,  A E  (also 
versus  A).  To  calculate  the  results,  a(t)  and  (3{t)  are  designed  as  if  A  =  0,  but  the  dynamics  is  carried  out  with  a 
non-zero  A,  in  particular  the  initial  state  is  the  actual  ground  state,  including  the  perturbation.  Note  that  A E 
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AE/(103/7W0) 


-4.4  -2.2  0  2.2  4.4 


A  (zN) 


Figure  4.  Excitation  energy  versus  different  tilt  values  of  the  external  potential  in  terms  of  the  energy  difference  between  both  wells 
(upper  axis)  and  values  of  the  A  parameter  (lower  axis),  when  using  the  1 1th  order  polynomial  in  the  evolution.  Same  parameters  as  in 
figure  2. 


should  be  more  than  a  thousand  vibrational  quanta  to  excite  the  final  energy  by  one  quantum.  In  [8]  an  energy 
increase  of  ten  phonons  at  about  150  zN  and  80  /is  separation  time  was  reported,  so  the  STA  ramps  definitely 
improve  the  robustness  against  bias. 

Further  experimental  limitations  maybe  due  to  random  fluctuations  in  the  potential  parameters,  or  higher 
order  terms  in  the  external  potential.  We  leave  these  important  issues  for  a  separate  study  but  note  that  the 
structure  of  the  STA  techniques  used  here  is  well  adapted  to  deal  with  noise  or  perturbations  [24-26]. 

Combining  STA  with  optimal- control- theory  methods  is  also  feasible,  see  e.g.  [12]. 

Finally,  we  compare  in  figure  5  the  performance  of  the  protocols  based  on  the  polynomials  (25)  and  (Al) 
with  a  simple  non- optimized  protocol  based  on  those  experimentally  used  in  [8].  There,  the  equilibrium 
distance  d  is  first  designed  as  d(t)  =  d(0)  +  [d(tf )  —  d(0)]s2  sin(s7r/2),  where  s  =  t/tf.  From  the  family  of 
possible  potential  ramps  consistent  with  this  function,  we  chose  a  polynomial  that  drives  a  from  a  (0)  =  a0  to 
a  ( tf )  =  —  a0/2  (as  in  figure  2)  and  whose  first  derivatives  are  0  at  both  boundary  times.  /3  is  given  by 
equation  (3).  For  the  times  analysed  in  figure  5,  the  method  based  on  equations  (25)  and  (Al)  clearly 
outperforms  the  non- optimized  ramp.  To  get  excitations  below  the  single  motional  excitation  quantum  with  the 
non-optimized  protocol,  final  times  as  long  as  tf  ~  80  /i s  would  be  needed,  which  is  in  line  with  current 
experiments. 

We  conclude  that  the  method  presented  here,  could  bring  a  clear  improvement  with  respect  to  the  best 
results  experimentally  reported  so  far  [7, 8].  The  parameters  required  are  realistic  in  current  trapped  ions 
laboratories.  The  simulations  show  that,  under  ideal  conditions,  the  separation  of  two  ions  could  be  performed 
in  a  few  oscillation  periods,  at  times  similar  to  those  required  for  other  operations  as  transport  [13]  or 
expansions  [14],  also  studied  with  STA. 
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Appendix.  Ansatz  for  p+ 

The  ansatz  for  p+  that  satisfies  the  BC  p+(0)  =  1,  p+(tf)  =  7+>  p+(%)  =  p+(t b)  =  p+(%)  ='p+(%)  =  0  with  two 
free  parameters  takes  the  form 

=  1  —  ^126  —  1267+  +  u\q  +  5an)s5 
T  ^420  —  42O7+  T  T  24 U\\^s^ 

—  ^540  —  540t+  T  IO^io  T  4 

+  (315  -  3157+  +  10aio  +  40an)s8 

—  ^70  —  7O7+  T  5aio  T"  15an)s9 

+  ai0s10  +  fln5n.  (Al) 

To  minimize  the  perturbation  energy  in  equation  (28),  three  free  parameters  are  introduced, 
p+  —  1  —  (l26  —  1267+  +  C10  +  5cn  +  15q2)s5 
+  ^420  —  42O7+  +  5qo  +  24cn  +  70ci2^s6 

-  (540  -  5407+  +  10cio  +  45cn  +  126ci2)s7 
T  ^315  —  31 57+  T  10cio  +  40cn  T  105q2^58 

—  ^70  —  7O7+  +  5cio  +  15cn  +  35q2^59 

+  Ci0510  +  CnS11  +  c!2s12.  (A2) 
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